Implementing the implicit Euler method for mass-spring systems
Mikko Kauppila aka uttumuttu of Faktory
1. Introduction
Mass-spring systems, mostly cloth simulators, have recently become popular effects in demos, thanks to the ever increasing processing capabilities. Indeed, decomposing a deformable system to point mass particles and springs is an intuitive, versatile and powerful model for simulating a great amount of physical phenomena. Numerous methods have been proposed for modeling the point masses and the springs, and for solving the cloth dynamics.
One of the groundbreaking works on this area is the work by Baraff and Witkin [Baraff98], Large Steps in Cloth Simulation, which proposes a stable integration method for solving the cloth dynamics, allowing to take vastly larger time steps in the simulation. The method is usually called the "backward Euler method" or the "implicit Euler method" (here we use the latter). While the method is elegant and powerful in itself, it's also considerably more difficult than explicit schemes. The aim of this paper is to present the implementation of the method instead of plain fancy math (efficiency is also given special attention). Math is, on the other hand, inevitable in the model, so a bit of familiarity with differential equations and linear algebra won't harm you.
1.1. Constructing a cloth model of points and particles
Although the implicit Euler method works for any kind of mass-spring systems, cloth (or more generally deformable surfaces) is probably the most popular of these. Therefore it seems natural to present how to construct the points and springs out of a three-dimensional polygonal cloth object. Here I shall follow the model of Baraff and Witkin [Baraff98] with slight changes.
The particles are simply the vertices of the polygonal object, each given an equal mass. An example: if the model constists of n vertices and the whole model is supposed to have a mass of m kilograms, then each particle is given a mass of m/n kilograms (this of course assumes that particles are at least approximately evenly distributed across the surface). This kind of equal mass assumption is not an absolute rule - it's just an easy and fairly accurate model.
For each edge in the object we construct a stretch spring between the particles forming the edge. The spring is given a stretch spring constant, a user-entered parameter which describes the stiffness of the spring. The greater the spring constant, the stiffer the spring.
We also need a model for the bending property of the cloth. For each edge (e) that connects two triangles (t1,t2), we construct a bending spring. The particles the spring connects are the vertices of the two triangles (t1,t2) that are not a part of the edge (e). This spring is associated with a bending spring constant, which should be far smaller than the stretch spring constant (for realistic look).
All springs (streching springs and bending springs) have a certain rest length which is a distance at which the particles of the spring want to stay from each other. It's determined for all springs in the beginning of the simulation by computing the original distance (the distance in the original, non-deformed cloth object) between the two particles connected by the spring.
2. The implicit Euler method - The formulas
In this section I'll present the formulas required for implementing the implicit Euler method. The cloth object has n particles. Each particle has a mass mi, a position xi, a velocity vi and a force fi associated with it (i being the index of the particle). Position, velocity and force of a particle are 3-vectors, while the mass is a scalar. We gather the positions, velocities and forces of all the particles into global 3n-vectors x, v, and f (where, for example, xi denotes the 3-vector presenting the position of the ith particle). The masses are gathered into a mass matrix M, which is a 3nx3n matrix, where each entry Mji (1 <= i,j <= n, or 0 <= ij <= n-1 for C programmers :) is a 3x3 matrix. The diagonal entries Mii are identity matrices multiplied by the corresponding particle's mass:
Mii = mi*I
All non-diagonal entries are zero matrices. Storing such matrix therefore only requires storing the diagonal part, which is memory-efficient. Denoting x,v as functions of time and f as a function of x and v it all boils down into the second-order differential equation:
x''(t) = M^-1 * f(x,v)
The mass matrix is just a fancy math way of expressing that the acceleration of a particle is the force acting on it divided by its mass. The inverse of the mass matrix is also easy to compute: the diagonal entries are simply inverted. The above equation can be separated in two first-order equations:
x'(t) = v(t)
v'(t) = M^-1 * f(x,v)
2.1. The implicit Euler method
A common problem in cloth simulation is that the simulator tends to "blow up". What's the cause of this instability? It's essentially caused by the spring forces, which are said to be stiff. The problem is that a normal Euler step, or any explicit method for simulating cloth dynamics, can propagate the forces only to a finite number of neighbour nodes in a single time step. An example of this is the following: consider a rope consisting of n particles, neighbouring particles connected with springs. When a strong pull is applied to one end of the rope, it will inevitably take n number of normal Euler steps before the opposite end of the string acknowledges the pull. In other words we say that the normal (explicit) Euler supports local interactions only. The idea behind the implicit Euler is that it forms a linear system of equations for computing the new state of the cloth, thus supporting global interactions naturally - and the opposite end of the string can feel our pull in a single time step.
There's a great introductionary text on implicit methods for differential equations available in [Baraff 2001]. The formulas here are exactly the same as in Baraff's paper [Baraff98]. The normal Euler method, usually called explicit Euler, discretizes the differential equation as:
dx = v * dt or alternatively dx = (v + dv)*dt dv = M^-1 * f(x0,v0)*dt
Where x0 and v0 are the positions and velocities at the beginning of the time step. This, however, is not stable for larger dt. The implicit Euler method uses the following discretization:
dx = (v + dv)*dt
dv = M^-1 * f(x0+dx,v0+dv)*dt
This essentially means that the force is computed at the end of the time step. Of course, the problem is that we do not know the positions and velocities at the end of the time step so we can't compute the force there either. More than that, the resulting system of equations in non-linear. Therefore we perform a step called linearization which approximates f(x0+dx,v0+dv) using the explicit Euler method:
f(x0+dx,v0+dv) = f0 + df/dx * dx + df/dv * dv
Here f0 = f(x0,v0). df/dx and df/dv are 3nx3n matrices called Jacobians. Their structure will be explained in section 2.2. When substituting this and solving for dv we get the following equation:
- dt*M^-1*df/dv - dt^2*M^-1*df/dx) * dv = dt*M^-1*(f0 + dt*df/dx*v0)
I is an identity matrix of size 3nx3n. As Baraff shows, this should be multiplied with M, otherwise the matrix on the left side won't be symmetric (if all particle masses aren't equal). We get a somewhat shorter formula:
- dt*df/dv - dt^2*df/dx) * dv = dt*(f0 + dt*df/dx*v0)
Compactly, what we're trying to solve is A*dv = b (where the matrix A and the vector b can be seen in the equation above). This is a linear system of equations - find a vector dv which multiplied by matrix A gives b as the result. In the implementation stage we solve this equation using the method of conjugate gradients (CG). dx is then trivial to compute given the vector dv.
2.2. Of forces and Jacobians
Our internal cloth dynamics consists of springs only. A spring (a stretch spring or a bending spring) is used to make two particle, i and j, "want" to be at a certain distance (the rest length) from each other. The force acting between the particles i and j is:
fij = (1-L/|xij|)*xij*k
Where xij = xi-xj, L is the rest length of the spring and k is the stiffness of the spring. This force is then added to fj and the opposite force (-fij) is added to fi. This is natural considering Newton's law of force and counterforce.
We can also use a model for damping the spring forces. Note that in implicit integration damping isn't absolutely necessary because the method will cause intrinsic damping to the solution. We define the damping force of a spring as:
dij = (vij)*kd
Where vij = vi-vj and kd is the damping constant (which should be quite small). This force is again added to fj and the opposite is added to fi.
So far we haven't discussed what the matrices df/dx and df/dv are. They are the special matrices called Jacobians. They are, again, 3nx3n matrices consisting of 3x3 submatrices. The submatrices are defined as follows:
(df/dx)ji = dfj/dxi
(df/dv)ji = dfj/dvi
These derivatives are nonzero only if particle i affects particle j. In other words, they are nonzero only if the particles i and j are connected with a spring. Therefore we define the Jacobian (Jxij for position, Jvij for velocity) of a spring constisting of particles ij as [Choi 2002]:
for position:
Jxij = { (xij*xij^t)/(xij^t*xij) + [I-(xij*xij^t)/(xij^t*xij)]*(1-L/|xij|) } * k
for velocity:
Jvij = I * kd
Where xij = xi-xj, L is the rest length of the spring and k is the stiffness of the spring. (xij^t*xij) denotes inner product, a simple dot product of xij with itself. (xij*xij^t) denotes outer product which results in a 3x3-matrix, defined as follows: let vector v=(a,b,c)^t. The outer product (v*v^t) is:
[ a*a a*b a*c ]
[ b*a b*b b*c ]
[ c*a c*b c*c ]
Jxij is added to dfj/dxi and dfi/dxj (the terms of the spring not on the diagonal) and subtracted from dfi/dxi and dfj/dxj (the terms of the spring on the diagonal). The same corresponding additions/subtractions are done for Jvij and df/dv entries.
As we shall see in the implementation section, we shall not need to store explicit matrix for storing these terms. Storing a one would be completely possible, though, because the matrix is sparse, that is, most of the entries are zero, so we only need to store the nonzero entries. In fact, because in practical cloth simulation the amount of springs that a particle is connected to is constant (around eight), the matrix requires only O(n) amount of storage. The matrix-free method was first proposed in [Volino 2000] - it's a great improvement to both the speed and the simplicity of the simulator.
3. The implicit integration - The implementation
In this section I'll outline all of the parts required for implementing a full implicit Euler method for mass-spring systems. First we'll cover basic structures (classes), and then the functions that operate with them.
A note on accuracy: it's recommended to use double precision floating point arithmetic for Vector3D and Matrix3x3 classes, and throughout the simulation. Especially the CG algorithm is heavily dependant on the floating point accuracy (more on this can be found in [Shewchuck 1996]).
3.1. Basic structures
Starting with a top-down approach, we could first create a Cloth class. It's a question of taste how the data in the class will be structured. The structures presented here are just examples.
class Cloth
{
public:
int nParticles; // number of particles
int nSprings; // number of springs
ClothParticle * particle; // the particle array (allocated dynamically)
ClothSpring * spring; // the spring array (allocated dynamically)
};
The dynamic allocation and loading of the particles, and generation of the springs, is not presented here, because it would take too much article space, and it depends completely of the implementation. Next we need to create the classes ClothParticle and ClothSpring.
class ClothParticle
{
public:
Vector3D x; // position
Vector3D v; // velocity
Vector3D F; // force
double m; // the mass
};
class ClothSpring
{
public:
int n[2]; // the indices of the particles the spring connects
Matrix3x3 Jx; // Jacobian with respect to position
Matrix3x3 Jv; // Jacobian with resepct to velocity
double r; // rest length
double ks; // stiffness constant (aka. spring constant)
double kd; // damping constant
};
3.2. Functions
3.2.1. Simulation, computing the forces and the Jacobians
The basic structures for the simulator aren't much of a joy before we write some functions that use them. Our main simulation function shall be void Cloth::Simulate( double dt ), which advances the cloth dynamics forward by time dt. Secondly we need a method that computes the forces acting on each particle (spring forces, including damping forces, gravity and external user-interaction forces), called void Cloth::ComputeForces(). We also define a function that computes the Jacobians for the springs. An example code for it is given here:
void Cloth::ComputeJacobians()
{
int i;
for( i=0; i<springs; i++ )
{
Vector3D dx = particle[spring[i].n[0]].x-particle[spring[i].n[1]].x;
Matrix3x3 dxtdx, I3x3;
dxtdx = OuterProduct( dx,dx );
I3x3.LoadIdentity();
double l = sqrt(dx*dx);
if( l != 0 ) l = 1.0/l;
dxtdx = dxtdx*(l*l);
spring[i].Jx = (dxtdx + (I3x3-dxtdx)*(1-spring[i].r*l))
* (spring[i].ks);
spring[i].Jv.LoadIdentity();
spring[i].Jv *= spring[i].kd;
}
3.2.2. Solving the linear system
Now that we have the basic code ready for computing forces and Jacobians, we can proceed to actually solving the new state of the cloth. It's accomplished by first computing the changes in velocities, dv, and then trivially computing the changes in positions, dx. The equation for solving dv was:
A * dv = b
where
A = M - dt*df/dv - dt^2*df/dx, a 3nx3n matrix
b = dt*(f0 + dt*df/dx*v0), a 3n vector
We begin by computing the b vector. The nontrivial task here is the multiplication df/dx*v0. For that one, we shall write a separate function called void Cloth::MultiplyDfDx( Vector3D * src, Vector3D * dst ). We exploit the special structure of df/dx, which is completely described by the Jacobians stored in the spring structures. The entries of df/dx a spring connecting particles i and j affects are ji,ij,ii and jj.
(df/dx)ji = (df/dx)ij = spring.Jx (df/dx)ii = (df/dx)jj = -spring.Jx
spring.Jx is naturally the Jacobian stored in the spring class. Therefore, the matrix multiplication dst=(df/dx)*src for spring connecting particles (i,j) is accomplished by doing the following additions/subtractions:
dst[i] += Jx*src[j];
dst[j] += Jx*src[i];
dst[i] -= Jx*src[i];
dst[j] -= Jx*src[j];
which even simplifies to:
dst[i] -= Jx*(src[i]-src[j]);
dst[j] += Jx*(src[i]-src[j]);
A code snippet:
void Cloth::MultiplyDfDx( Vector3D * src, Vector3D * dst )
{
int i;
for( i=0; i<particles; i++ )
dst[i] = Vector3D(0,0,0); // initialization
for( i=0; i<springs; i++ )
{
Vector3D temp = spring[i].Jx*(src[spring[i].n[0]]
-src[spring[i].n[1]]);
dst[spring[i].n[0]] -= temp;
dst[spring[i].n[1]] += temp;
}
}
Next we need a procedure that multiplies any vector (Vector3D src[particles]) with the matrix A. This is required for the conjugate gradient (CG) algorithm that solves the final equation. The procedure consists of three parts (multiplying by M, multiplying by -dt*df/dv and finally multiplying by -dt^2*df/dx, and finally summing these up). The codes for multiplying by df/dv and df/dx are exactly similar to the void Cloth::MultiplyDfDx( Vector3D * src, Vector3D * dst ), so we won't repeat the code here.
After this comes the CG algorithm. The method of conjugate gradients solves a linear system of equations (A * dv = b) iteratively, that is, by first making an initial guess and then correcting this guess on each iteration. The algorithm works best for symmetric and positive-definite matrices (such as A, but see section 4.4 for further details on its positive-definitiness), and only requires multiplications through A (which makes it ideal when A is sparse). A great introduction and deep analysis of the method can be found on [Shewchuck 1996]. Here's the algorithm in pseudo-code (following Shewchuck's paper):
dv <= 0 // null vector is a good choice for initial guess
i <= 0 // a scalar
r <= b - A*dv // a vector
d <= r
epsNew <= r^t*r // dot product, a scalar
eps0 = epsNew
while( i < iMax && epsNew > eps*eps0 )
{
q <= Ad
alpha = epsNew / (d^t*q)
dv <= dv + alpha*d
r <= r - alpha*q
epsOld = epsNew
epsNew = r^t*r
beta <= epsNew/epsOld
d <= r + beta*d
i <= i+1
}
eps is a user-defined parameter (eps < 1) that defines when the answer (dv) is accurate enough. The smaller the value of eps, the more accurate the answer, but the larger the number of iterations taken. A and b are, of course, the matrix and vector of the A*dv=b problem.
3.2.3 Adaptive time stepping
Even though the implicit integration is inherently stable, there are still times when the time step must be decreased. This is especially true in collision with large energies. As [Bridson 2002], we use the "rule of thumb" of computational mechanics: after a time step has been taken, we check if any spring has deformed more than 10 percents (15-20 percents goes well in practical simulation, too). If it has, then we discard the new state, restore the old state and halve the time step. After two subsequent succesful steps we double the time step again. This is a bit similar to the adaptive time stepping of [Baraff 1998], although Baraff's one is somewhat more complex in order to reduce the number of needless tests for larger time steps.
4. Further goodies
In this section further optimizations and enhancements to the simulator are presented. The subject is large, and everything simply can't be discussed here, but I'll try to give a brief review of some methods.
4.1. Optimizing the matrix-vector multiplication
Most of the time in the simulation is spent in the conjugate gradient algorithm. The most time consuming part of the CG algorithm itself is the matrix-vector multiplication by A. Therefore we should seek for strategies for minimizing the cost of this matrix-vector multiplication.
A = M - dt*df/dv - dt^2*df/dx
We split A in two parts:
A = M+B
B = -dt*df/dv - dt^2*df/dx
The submatrices of B can be precomputed for each spring, just like the Jacobian entries are computed. Notice, that the 3x3-submatrices are symmetric, so we don't need to store all the nine values, but only six of them. Storing a lot of Jacobian information in a spring structure kills the cache performance. A good speedup is achieved by defining a precomputed, cache-friendly class for the springs that is used in void Cloth::MultiplyMatrix( Vector3D * src, Vector3D * dst ). An example would be:
class ClothSpringPrecomputed
{
public:
int n[2]; // the indices of the particles connected by the spring
double m[6]; // the precomputed Jacobian (-dt*Jv - dt^2*Jx)
// because of symmetricity of the Jacobians, only six values
// need to be stored
};
4.2. Preconditioning the conjugate gradient algorithm
The conjugate gradient algorithm iteratively finds an approximate solution for a linear system of equations, in our case A*dv=b. The number of iterations taken to produce a reasonably accurate solution can be decreased through a technique called preconditioning. There's a general consensus that the Jacobi preconditioner is a great and simple choice for cloth simulation [Shewchuck 1996, Baraff 1998, Choi 2002]. The interested reader is pointed to those references.
4.3. Collision handling
Various types of collision detection/handling algorithms have been proposed in the literature. Baraff's original work [Baraff 1998] used a clever constraint mechanism with the CG algorithm (for cloth/solid collision and contact) and penalty forces (for cloth/cloth collision and contact). A completely different kind of collision handling algorithm can be found in the name monster paper Robust Treatment of Collision, Contact and Friction for Cloth Modeling [Bridson, Fedkiw and Anderson 2002].
5. Acknowledgements and conclusions
Special thanks go to Kwang-Jin Choi for his help via email and to Safe Zhou for his careful proofreading.
In this tutorial, an implementation of the implicit Euler method for mass-spring systems was presented. I hope the tutorial is a help to those of us (like me) who found the original paper of Baraff somewhat confusing. Feedback and questions are warmly welcomed at mijoka02@oppi.edu.ouka.fi.
References
Baraff, D. & Witkin, A. Large Steps in Cloth Simulation. In Proceedings of SIGGRAPH, 1998
Baraff, D. Implicit Methods for Differential Equations. SIGGRAPH course notes, 2001
Bridson, R., Fedkiw, R. and Anderson, J. Robust Treatment of Collision, Contact and Friction for Cloth Animation. In Proceedings of SIGGRAPH, 2002
Choi, K.-J. & Ko, H.-S. Stable but Responsive Cloth. In Proceedings of SIGGRAPH, 2002
Volino, P. & Magnenat-Thalmann, N. Implementing Fast Cloth Simulation with Collision Response. Computer Graphics International 2000